import numpy as np
import pandas as pd
from osgeo import gdal


# 根据像元行列号计算经纬度
def pixel2coord(geo_transform, col, row):
    x_origin, pixel_width, x_rotation, y_origin, y_rotation, pixel_height = geo_transform
    lon = x_origin + pixel_width * col
    lat = y_origin + pixel_height * row
    return lon, lat


# 根据经纬度计算像元行列号
def coord2pixel(geo_transform, lon, lat):
    col = int((lon - geo_transform[0]) / geo_transform[1])
    row = int((lat - geo_transform[3]) / geo_transform[5])
    return col, row

def read2PSM(station_path, PSM_path, sm_dir):
    station_df = pd.read_csv(station_path)
    psm_df = pd.read_excel(PSM_path, sheet_name="水分均值")

    # 根据日期计算儒略日
    # targetday = datetime.date(int(str(date)[0:4]), int(str(date)[4:6]), int(str(date)[6:8]))
    # day = targetday - datetime.date(targetday.year - 1, 12, 31)
    # julian_day = str(date)[0:4] + str(day.days)

    for i, psm in psm_df.iterrows():
        station_id = psm['station']
        # print(psm['date'].strftime("%Y%m%d"))
        date = psm['date'].strftime("%Y%m%d")
        # print(station_id)

        sm_value_list = []
        for j, station in station_df.iterrows():
            if station['station'].split("-")[0] == str(station_id):
                # print(station['station'])
                lon = station['longitude']
                lat = station['latitude']
                sm_path = sm_dir + "\\" + date + ".tif"
                sm_dataset = gdal.Open(sm_path)
                sm_geo = sm_dataset.GetGeoTransform()
                sm_data = sm_dataset.ReadAsArray()
                col, row = coord2pixel(sm_geo, lon, lat)
                sm_value = sm_data[row, col]
                # print(sm_value)
                if 0 <= sm_value <= 1:
                    sm_value_list.append(sm_value)
        psm_df.loc[i, 'date'] = str(date)
        if len(sm_value_list) > 0:
            print(sm_value_list)
            psm_df.loc[i, 'mcca'] = np.mean(sm_value_list)
    print(psm_df)
    psm_df.to_excel(r"G:\test\SMAP_NEW\PSM_RCD_SM_mean.xlsx", index=False)


if __name__ == '__main__':
    station_path = r"G:\test\SMAP_NEW\station_id.csv"
    psm_path = r"G:\test\SMAP_NEW\PSM_RCD_SM1.xlsx"
    sm_dir = r"G:\test\SMAP_NEW\MCCA\TIF\2016"
    read2PSM(station_path, psm_path, sm_dir)

